rm(list = ls())
require(zTree)
require(stargazer)
require(mfx)
require(stats4)
require(nlWaldTest)
require("multiwayvcov")
require(msm)
require(miceadds)


zTT1= zTreeTables("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/Ambueletalexperiments/Data/August2Expansionpath2021.xls")
subjects1=zTT1$subjects

zTT2= zTreeTables("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/Ambueletalexperiments/Data/September21ExpansionPath2021.xls")
subjects2=zTT2$subjects
subjects2$Subject=subjects2$Subject+100

zTT3= zTreeTables("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/Ambueletalexperiments/Data/September 29ExpansionPath2021.xls")
subjects3=zTT3$subjects
subjects3$Subject=subjects3$Subject+200

zTT4= zTreeTables("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/Ambueletalexperiments/Data/October13expansionpath2021.xls")
subjects4=zTT4$subjects
subjects4$Subject=subjects4$Subject+300


subjects=rbind(subjects1,subjects2,subjects3,subjects4)

#Clean data
subjects$responder=(subjects$"response1[50]">0)
subjects=subset(subjects,responder>0)

uidlist=unique(subjects$Subject)
incentivec=c(5,95)

#create a frame with observations
obsframe=data.frame(uid=vector('numeric'),block=vector('numeric'), incentive=vector('numeric'), number=vector('numeric'), state=vector('numeric'), response=vector('numeric'),react=vector('numeric'))
for(i in 1:length(uidlist)){
playframe=subjects[i,]
uidplace=uidlist[i]

#block 1
for (j in 1:50){
blockplace=1
incentiveplace=incentivec[playframe[1,683]]
numberplace=j
stateplace=playframe[1,(736+j)]
responseplace=playframe[1,(72+j)]
reactplace=playframe[1,(122+j)]
addframe=data.frame(uid=uidplace,block=blockplace, incentive=incentiveplace, number=numberplace, state=stateplace, response=responseplace,react=reactplace)
obsframe=rbind(obsframe,addframe)
}

#block 2
for (j in 1:50){i
blockplace=2
incentiveplace=incentivec[playframe[1,684]]
numberplace=j
stateplace=playframe[1,(786+j)]
responseplace=playframe[1,(222+j)]
reactplace=playframe[1,(272+j)]
addframe=data.frame(uid=uidplace,block=blockplace, incentive=incentiveplace, number=numberplace, state=stateplace, response=responseplace,react=reactplace)
obsframe=rbind(obsframe,addframe)

}

}



#define correctness
obsframe$correct=(obsframe$state==obsframe$response)

#mean correctness
correct5=sum(obsframe$correct*(obsframe$incentive==5))/sum(obsframe$incentive==5)
correct95=sum(obsframe$correct*(obsframe$incentive==95))/sum(obsframe$incentive==95)
correctness=c(correct5,correct95)

obsframe5=subset(obsframe,incentive==5)
obsframe95=subset(obsframe,incentive==95)

#change terminology to work with older scripts
expframe=obsframe

expframe$incentive5=(expframe$incentive==5)
expframe$incentive95=(expframe$incentive==95)

expframe$Bresponse <-expframe$response==2
expframe$Bstate <- expframe$state==2
expframe$Aresponse <-expframe$response==1
expframe$Astate <- expframe$state==1

expframe5=subset(expframe,incentive==5)
expframe95=subset(expframe,incentive==95)

#NIAS in aggregate

#NIAS in aggregate
#Base Regression
expframe$incentivedec=expframe$incentive/100
modelNIAS=lm(Aresponse~factor(Astate+incentivedec)+0,data=expframe)

#Adjust COV
vcovmodNIAScluster=cluster.vcov(modelNIAS, expframe$uid)

modelniassterr=vector("numeric")
#Model Nias sterrors
for(i in 1:4){
modelniassterr[i]=sqrt(vcovmodNIAScluster[i,i])
}


Textstr=vector("character")
Pvalvec=vector("numeric")

#Tests
Textstr[1]="(b[3]-b[1])"
wald1=nlWaldtest(modelNIAS,texts=Textstr,Vcov=vcovmodNIAScluster)
Pvalvec[1]=wald1$p.value

Textstr[1]="(b[4]-b[2])"
wald2=nlWaldtest(modelNIAS,texts=Textstr,Vcov=vcovmodNIAScluster)
Pvalvec[2]=wald1$p.value

Textstr=vector("character")
#Joint test
Textstr[1]="(b[3]-b[1])"
Textstr[2]="(b[4]-b[2])"
nlWaldtest(modelNIAS,texts=Textstr,Vcov=vcovmodNIAScluster)

pAAvec=c(modelNIAS$coefficient[3],modelNIAS$coefficient[4])
pABvec=c(modelNIAS$coefficient[1],modelNIAS$coefficient[2])



#Incentive level 5
expframe5=subset(expframe,incentive==5)
model5=logitmfx(Aresponse~Astate,data=expframe5,robust=TRUE, clustervar1 = "uid")
prob5=model5$mfxest[1,4]
state1achance5=sum(expframe5$Aresponse*expframe5$Astate)/sum(expframe5$Astate)
state2achance5=sum(expframe5$Aresponse*expframe5$Bstate)/sum(expframe5$Bstate)

#Incentive level 95
expframe95=subset(expframe,incentive==95)
model95=logitmfx(Aresponse~Astate,data=expframe95,robust=TRUE, clustervar1 = "uid") 
prob95=model95$mfxest[1,4]
state1achance95=sum(expframe95$Aresponse*expframe95$Astate)/sum(expframe95$Astate)
state2achance95=sum(expframe95$Aresponse*expframe95$Bstate)/sum(expframe95$Bstate)


ANIASreport=data.frame(Incentive_Level=incentivec, P_Aresp_in_Astate=pAAvec,P_Aresp_in_Bstate=pABvec,Pval=Pvalvec)


#***********************************************
stargazer(ANIASreport,summary=FALSE,rownames=FALSE)
#***********************************************

